library("dplyr")
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library("lubridate")
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library("MASS")
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
library("MuMIn")
library("forecast")
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
source("plot.roots.R")
data <- read.csv("unemployment.csv", header = TRUE)
data <- data[data$State.Area == "California", ]

data <- data %>%
  mutate(Total.Unemployment.in.State.Area = as.numeric(gsub("[^0-9]", "", Total.Unemployment.in.State.Area))) 

#Subset training to till 2018, testing till 2019
data_training <- data %>% 
  filter(Year >= 2010 & Year <= 2018)

data_testing <- data %>% 
  filter(Year >= 2010 & Year <= 2019)

#Make time series
data_training_ts = ts(data_training$Total.Unemployment.in.State.Area, start = c(2010, 1), end = c(2018, 1), frequency = 12) 

data_testing_ts = ts(data_testing$Total.Unemployment.in.State.Area, start = c(2010, 1), end = c(2019, 1), frequency = 12)

#Plot time series
plot(data_training_ts, xlab = "Year", ylab = "Unemployment", main = "Figure 1: Unemployment in California from 2010 to 2018", ylim = c(min(data_training_ts), max(data_training_ts)))

#Let's transform data to make it more normal

#Box-Cox Transformation
bcTransform <- boxcox(data_training_ts ~ as.numeric(1:length(data_training_ts)))

bcTransform$x[which(bcTransform$y == max(bcTransform$y))]
## [1] 0.3838384
lambda = bcTransform$x[which(bcTransform$y == max(bcTransform$y))]
data_ts.bc = (1/lambda)*(data_training_ts^lambda-1)
plot(data_ts.bc, main = "Box-Cox Transformation of Unemployment")

hist(data_ts.bc, col = "lightblue", xlab = "Year", ylab = "Unemployment", main = "Box-Cox Transformation Histogram")

#Log Transformation
data_ts.log <- log(data_training_ts)
plot(data_ts.log, main = "Log Transformation of Unemployment")

hist(data_ts.log, col="lightblue", xlab = "Year", ylab = "Unemployment", main = "Log Transformation Histogram")

#Sqrt Transformation
data_ts.sqrt <- sqrt(data_training_ts)
plot(data_ts.sqrt, main = "SQRT Transformation of Unemployment")

hist(data_ts.sqrt, col="lightblue", xlab = "Year", ylab = "Unemployment", main = "Figure 2: SQRT Transformation Histogram")

#The Box-Cox command gives us lambda = 0.3838384. Lambda = 0 corresponds to a Log transformation, while lambda = 0.5 corresponds to a SQRT transformation. Since 0.5 is closer to 0.3838384 and is within the 95% confidence interval, we choose the SQRT transformation. This is also supported by a more symmetric looking histogram.
#Difference to get rid of trend (lag 1)
diff_lag_1_data_ts.transformed = diff(data_ts.sqrt, 1)
plot(diff_lag_1_data_ts.transformed, xlab = "Year", ylab = "Unemployment", main ="De-Trended Unemployment")
abline(h = mean(diff_lag_1_data_ts.transformed),lty = 2)

var(diff_lag_1_data_ts.transformed)
## [1] 13.35966
#Difference to get rid of seasonality (lag 1 + lag 12)
diff_lag_1_and_12_data_ts.transformed = diff(diff_lag_1_data_ts.transformed, 12)
plot(diff_lag_1_and_12_data_ts.transformed, xlab = "Year", ylab = "Unemployment", main ="De-Trended and De-seasonalized Unemployment")
abline(h = mean(diff_lag_1_and_12_data_ts.transformed),lty = 2)

var(diff_lag_1_and_12_data_ts.transformed)
## [1] 23.9294
#Differencing at lag 1 gives us a variance of 3.822456e-05, while differencing at lag 1 and lag 12 gives us a variance of 6.681259e-05. Since the former is a lower variance, we choose to de-trend only.
#De-Trended ACF and PACF

#ACF
acf(diff_lag_1_data_ts.transformed, xlim = c(0, 5), main = " Figure 3: ACF of Unemployment in California from 2010 to 2018", col = "orchid", xlab = "Lag")

#PACF
pacf(diff_lag_1_data_ts.transformed, xlim = c(0, 5), main = "Figure 4: PACF of Unemployment in California from 2010 to 2018", col = "orchid", xlab = "Lag")

#ARIMA looks like the most appropriate model. p is most likely either 2 or 4 since there is a strong AR component, d is either 1 or 2 since we differenced at 1 to remove trend but might need to account for more differencing, and q is definitely 0 since there is no MA component (seasonality).
#ARIMA testing
arima(data_ts.sqrt, order=c(2,1,0), method="ML")
## 
## Call:
## arima(x = data_ts.sqrt, order = c(2, 1, 0), method = "ML")
## 
## Coefficients:
##          ar1      ar2
##       1.6937  -0.7366
## s.e.  0.0687   0.0690
## 
## sigma^2 estimated as 1.144:  log likelihood = -144.97,  aic = 295.93
arima(data_ts.sqrt, order=c(4,1,0), method="ML")
## 
## Call:
## arima(x = data_ts.sqrt, order = c(4, 1, 0), method = "ML")
## 
## Coefficients:
##          ar1      ar2     ar3      ar4
##       2.3097  -2.1781  0.9204  -0.0640
## s.e.  0.1055   0.2472  0.2479   0.1063
## 
## sigma^2 estimated as 0.4523:  log likelihood = -101.74,  aic = 213.48
arima(data_ts.sqrt, order=c(2,2,0), method="ML")
## 
## Call:
## arima(x = data_ts.sqrt, order = c(2, 2, 0), method = "ML")
## 
## Coefficients:
##          ar1      ar2
##       1.2784  -0.7952
## s.e.  0.0619   0.0613
## 
## sigma^2 estimated as 0.4604:  log likelihood = -99.31,  aic = 204.62
arima(data_ts.sqrt, order=c(4,2,0), method="ML")
## 
## Call:
## arima(x = data_ts.sqrt, order = c(4, 2, 0), method = "ML")
## 
## Coefficients:
##          ar1      ar2     ar3      ar4
##       1.3482  -1.1001  0.3960  -0.2519
## s.e.  0.1040   0.1747  0.1759   0.1046
## 
## sigma^2 estimated as 0.4315:  log likelihood = -96.39,  aic = 202.78
AICc(arima(data_ts.sqrt, order=c(2,1,0), method="ML"))
## [1] 296.1931
AICc(arima(data_ts.sqrt, order=c(4,1,0), method="ML"))
## [1] 214.1453
AICc(arima(data_ts.sqrt, order=c(2,2,0), method="ML"))
## [1] 204.8874
AICc(arima(data_ts.sqrt, order=c(4,2,0), method="ML"))
## [1] 203.4565
#The AICc values for the models are as follows: ARIMA(2,1,0) = 296.1931, ARIMA(4,1,0) = 214.1453, ARIMA(2,2,0) = 204.8874, and ARIMA(4,2,0) = 203.4565. Since ARIMA(4,2,0) has the lowest AICc value, we select this model.
#Let's check for stationarity since we know it is invertible (definition of AR)
#ARIMA(2,1,0)
coeff1 <- c(1, -1.6937, 0.7366)
roots1 <- polyroot(coeff1)
plot.roots(NULL, roots1, main = "Roots of AR Polynomial for ARIMA(2,1,0)")

#roots are outside the unit circle, so our model is stationary

#ARIMA(4,1,0)
coeff2 <- c(1, -2.3097, 2.1781, -0.9204, 0.0640)
roots2 <- polyroot(coeff2)
plot.roots(NULL, roots2, main = "Roots of AR Polynomial for ARIMA(4,1,0)")

#roots are outside the unit circle, so our model is stationary

#ARIMA(2,2,0)
coeff3 <- c(1, -1.2784, 0.7952)
roots3 <- polyroot(coeff3)
plot.roots(NULL, roots3, main = "Roots of AR Polynomial for ARIMA(2,2,0)")

#roots are outside the unit circle, so our model is stationary

#ARIMA(4,2,0)
coeff4 <- c(1, -1.3482, 1.1001, -0.3960, 0.2519)
roots4 <- polyroot(coeff4)
plot.roots(NULL, roots4, main = "Roots of AR Polynomial for ARIMA(4,2,0)")

#roots are outside the unit circle, so our model is stationary
#Let's fit the model ARIMA(2,1,0) and get residuals
fit <- arima(data_ts.sqrt, order=c(2,1,0), method="ML")
res1 <- residuals(fit)

#Let's look at the residuals closely
hist(res1, density=20, breaks=20, col="pink", prob=TRUE)
mean <- mean(res1)
standard_deviation <- sqrt(var(res1))
curve(dnorm(x, mean, standard_deviation), add=TRUE )

plot.ts(res1, main = "ARIMA(2,1,0) Residuals")

#Q-Q Plot
fitt <- lm(res1 ~ as.numeric(1:length(res1))); abline(fitt, col="purple")
abline(h=mean(res1), col="lightblue")

qqnorm(res1,main= "Q-Q Plot")
qqline(res1,col="darkblue")

#Seems to follow normal distribution.

#ACF and PACF
acf(res1, lag.max=40)

pacf(res1, lag.max=40)

acf(res1^2, lag.max=40)

#Shapiro-Wilk Test
shapiro.test(res1)
## 
##  Shapiro-Wilk normality test
## 
## data:  res1
## W = 0.98112, p-value = 0.1771
#Since p-value = 0.1771, we fail to reject the null hypothesis that the residuals are normally distributed. Therefore the residuals are normally distributed.

#Box-Pierce Test
Box.test(res1, lag = 12, type = c("Box-Pierce"), fitdf = 2)
## 
##  Box-Pierce test
## 
## data:  res1
## X-squared = 113.93, df = 10, p-value < 2.2e-16
#Since p-value < 2.2e-16, we reject the null hypothesis that the residuals are white noise. Therefore the residuals are not white noise and have some autocorrelation.

#Ljung-Box Test
Box.test(res1, lag = 12, type = c("Ljung-Box"), fitdf = 2)
## 
##  Box-Ljung test
## 
## data:  res1
## X-squared = 122.6, df = 10, p-value < 2.2e-16
#Since p-value < 2.2e-16, we reject the null hypothesis that the residuals are independently distributed. Therefore the residuals are not independent and there is significant autocorrelation.

#McLeod-Li Test
Box.test(res1^2, lag = 12, type = c("Ljung-Box"), fitdf = 0)
## 
##  Box-Ljung test
## 
## data:  res1^2
## X-squared = 11.826, df = 12, p-value = 0.4597
#Since p-value = 0.4597, we fail to reject the null hypothesis that the variance of the residuals is white noise. Therefore there is no significant autocorrelation in the variance.

#Yule-Walker Test
ar(res1, aic = TRUE, order.max = NULL, method = c("yule-walker"))
## 
## Call:
## ar(x = res1, aic = TRUE, order.max = NULL, method = c("yule-walker"))
## 
## Coefficients:
##       1        2        3        4        5        6        7        8  
##  0.5198  -0.4311  -0.1273  -0.3891   0.0915  -0.1138  -0.0702  -0.0515  
##       9       10       11       12       13  
##  0.2757   0.0337   0.0415  -0.0844   0.2884  
## 
## Order selected 13  sigma^2 estimated as  0.4928
#sigma^2 = 0.4928.

#This model fails the Box-Pierce Test and Ljung-Box Test.
#Let's fit the model ARIMA(4,1,0) and get residuals
fit <- arima(data_ts.sqrt, order=c(4,1,0), method="ML")
res2 <- residuals(fit)

#Let's look at the residuals closely
hist(res2, density=20, breaks=20, col="pink", prob=TRUE)
mean <- mean(res2)
standard_deviation <- sqrt(var(res2))
curve(dnorm(x, mean, standard_deviation), add=TRUE )

plot.ts(res2, main = "ARIMA(4,1,0) Residuals")

#Q-Q Plot
fitt <- lm(res2 ~ as.numeric(1:length(res2))); abline(fitt, col="purple")
abline(h=mean(res2), col="lightblue")

qqnorm(res2,main= "Q-Q Plot")
qqline(res2,col="darkblue")

#Seems to follow normal distribution.

#ACF and PACF
acf(res2, lag.max=40)

pacf(res2, lag.max=40)

acf(res2^2, lag.max=40)

#Shapiro-Wilk Test
shapiro.test(res2)
## 
##  Shapiro-Wilk normality test
## 
## data:  res2
## W = 0.9904, p-value = 0.7159
#Since p-value = 0.7159, we fail to reject the null hypothesis that the residuals are normally distributed. Therefore the residuals are normally distributed.

#Box-Pierce Test
Box.test(res2, lag = 12, type = c("Box-Pierce"), fitdf = 2)
## 
##  Box-Pierce test
## 
## data:  res2
## X-squared = 25.517, df = 10, p-value = 0.004448
#Since p-value = 0.004448, we reject the null hypothesis that the residuals are white noise. Therefore the residuals are not white noise and have autocorrelation.

#Ljung-Box Test
Box.test(res2, lag = 12, type = c("Ljung-Box"), fitdf = 2)
## 
##  Box-Ljung test
## 
## data:  res2
## X-squared = 28.573, df = 10, p-value = 0.00146
#Since p-value = 0.00146, we reject the null hypothesis that the residuals are independently distributed. Therefore the residuals are not independent and there is significant autocorrelation.

#McLeod-Li Test
Box.test(res2^2, lag = 12, type = c("Ljung-Box"), fitdf = 0)
## 
##  Box-Ljung test
## 
## data:  res2^2
## X-squared = 19.714, df = 12, p-value = 0.0727
#Since p-value = 0.0727, we fail to reject the null hypothesis that the variance of the residuals is white noise. Therefore there is no significant autocorrelation in the variance.

#Yule-Walker Test
ar(res2, aic = TRUE, order.max = NULL, method = c("yule-walker"))
## 
## Call:
## ar(x = res2, aic = TRUE, order.max = NULL, method = c("yule-walker"))
## 
## 
## Order selected 0  sigma^2 estimated as  0.4669
#sigma^2 = 0.4669.

#This model fails the Box-Pierce Test and Ljung-Box Test.
#Let's fit the model ARIMA (2,2,0) and get residuals
fit <- arima(data_ts.sqrt, order=c(2,2,0), method="ML")
res3 <- residuals(fit)

#Let's look at the residuals closely
hist(res3, density=20, breaks=20, col="pink", prob=TRUE)
mean <- mean(res3)
standard_deviation <- sqrt(var(res3))
curve(dnorm(x, mean, standard_deviation), add=TRUE )

plot.ts(res3, main = "ARIMA(2,2,0) Residuals")

#Q-Q Plot
fitt <- lm(res3 ~ as.numeric(1:length(res3))); abline(fitt, col="purple")
abline(h=mean(res3), col="lightblue")

qqnorm(res3,main= "Q-Q Plot")
qqline(res3,col="darkblue")

#Seems to follow normal distribution.

#ACF and PACF
acf(res3, lag.max=40)

pacf(res3, lag.max=40)

acf(res3^2, lag.max=40)

#Shapiro-Wilk Test
shapiro.test(res3)
## 
##  Shapiro-Wilk normality test
## 
## data:  res3
## W = 0.99178, p-value = 0.82
#Since p-value = 0.82, we fail to reject the null hypothesis that the residuals are normally distributed. Therefore the residuals are normally distributed.

#Box-Pierce Test
Box.test(res3, lag = 12, type = c("Box-Pierce"), fitdf = 2)
## 
##  Box-Pierce test
## 
## data:  res3
## X-squared = 22.382, df = 10, p-value = 0.01327
#Since p-value = 0.01327 we reject the null hypothesis that the residuals are white noise. Therefore the residuals are not white noise and have some autocorrelation.

#Ljung-Box Test
Box.test(res3, lag = 12, type = c("Ljung-Box"), fitdf = 2)
## 
##  Box-Ljung test
## 
## data:  res3
## X-squared = 25.184, df = 10, p-value = 0.005008
#Since p-value = 0.005008, we reject the null hypothesis that the residuals are independently distributed. Therefore the residuals are not independent and there is significant autocorrelation.

#McLeod-Li Test
Box.test(res3^2, lag = 12, type = c("Ljung-Box"), fitdf = 0)
## 
##  Box-Ljung test
## 
## data:  res3^2
## X-squared = 13.144, df = 12, p-value = 0.3587
#Since p-value = 0.3587, we fail to reject the null hypothesis that the variance of the residuals is white noise. Therefore there is no significant autocorrelation in the variance.

#Yule-Walker Test
ar(res3, aic = TRUE, order.max = NULL, method = c("yule-walker"))
## 
## Call:
## ar(x = res3, aic = TRUE, order.max = NULL, method = c("yule-walker"))
## 
## 
## Order selected 0  sigma^2 estimated as  0.4993
#sigma^2 = 0.4993

#This model fails the Box-Pierce Test and Ljung-Box Test.
#Let's fit the model ARIMA(4,2,0) and get residuals
fit <- arima(data_ts.sqrt, order=c(4,2,0), method="ML")
res4 <- residuals(fit)

#Let's look at the residuals closely
hist(res4, density=20, breaks=20, col="pink", prob=TRUE)
mean <- mean(res4)
standard_deviation <- sqrt(var(res4))
curve(dnorm(x, mean, standard_deviation), add=TRUE )

plot.ts(res4, main = "ARIMA(4,2,0) Residuals")

#Q-Q Plot
fitt <- lm(res4 ~ as.numeric(1:length(res4))); abline(fitt, col="purple")
abline(h=mean(res4), col="lightblue")

qqnorm(res4,main= "Q-Q Plot")
qqline(res4,col="darkblue")

#Seems to follow normal distribution.

#ACF and PACF
acf(res4, lag.max=40)

pacf(res4, lag.max=40)

acf(res4^2, lag.max=40)

#Shapiro-Wilk Test
shapiro.test(res4)
## 
##  Shapiro-Wilk normality test
## 
## data:  res4
## W = 0.98406, p-value = 0.2903
#Since p-value = 0.2903, we fail to reject the null hypothesis that the residuals are normally distributed. Therefore the residuals are normally distributed.

#Box-Pierce Test
Box.test(res4, lag = 12, type = c("Box-Pierce"), fitdf = 2)
## 
##  Box-Pierce test
## 
## data:  res4
## X-squared = 11.622, df = 10, p-value = 0.3111
#Since p-value = 0.3111, we fail to reject the null hypothesis that the residuals are white noise. Therefore the residuals are white noise and have no autocorrelation.

#Ljung-Box Test
Box.test(res4, lag = 12, type = c("Ljung-Box"), fitdf = 2)
## 
##  Box-Ljung test
## 
## data:  res4
## X-squared = 13.133, df = 10, p-value = 0.2163
#Since p-value = 0.2163, we fail to reject the null hypothesis that the residuals are independently distributed. Therefore the residuals are independent and there is no significant autocorrelation.

#McLeod-Li Test
Box.test(res4^2, lag = 12, type = c("Ljung-Box"), fitdf = 0)
## 
##  Box-Ljung test
## 
## data:  res4^2
## X-squared = 10.735, df = 12, p-value = 0.5518
#Since p-value = 0.5518, we fail to reject the null hypothesis that the variance of the residuals is white noise. Therefore there is no significant autocorrelation in the variance.

#Yule-Walker Test
ar(res4, aic = TRUE, order.max = NULL, method = c("yule-walker"))
## 
## Call:
## ar(x = res4, aic = TRUE, order.max = NULL, method = c("yule-walker"))
## 
## 
## Order selected 0  sigma^2 estimated as  0.4698
#sigma^2 = 0.4698

#This model passes all tests.
#Since ARIMA(4,2,0) is the only model that checks all the diagnostic criteria and has the lowest AICc value, we are confident this is the appropriate model. 

#Equation for selected model ARIMA(4,2,0)
model <- arima(data_ts.sqrt, order=c(4,2,0), method="ML")
summary(model)
## 
## Call:
## arima(x = data_ts.sqrt, order = c(4, 2, 0), method = "ML")
## 
## Coefficients:
##          ar1      ar2     ar3      ar4
##       1.3482  -1.1001  0.3960  -0.2519
## s.e.  0.1040   0.1747  0.1759   0.1046
## 
## sigma^2 estimated as 0.4315:  log likelihood = -96.39,  aic = 202.78
## 
## Training set error measures:
##                       ME      RMSE       MAE          MPE       MAPE       MASE
## Training set -0.06814014 0.6852723 0.5639344 -0.003994467 0.04696552 0.08781038
##                     ACF1
## Training set -0.01794226
#Δ_2 sqrt(U_t) = (1 - 1.3482B + 1.1001B^2 - 0.3960B^3 + 0.2519B^4)Z_t
#Sqrt transformed data plot with 12 forecasts
fit <- arima(sqrt(data_training_ts), order=c(4,2,0), method="ML")
forecast(fit)
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Feb 2018       913.6224 912.7806 914.4642 912.3349 914.9098
## Mar 2018       909.6700 906.7285 912.6115 905.1713 914.1686
## Apr 2018       905.8627 899.7143 912.0112 896.4596 915.2659
## May 2018       901.3192 891.3253 911.3132 886.0348 916.6036
## Jun 2018       895.6559 881.6620 909.6498 874.2540 917.0577
## Jul 2018       889.0394 871.2290 906.8499 861.8008 916.2781
## Aug 2018       882.0419 860.7158 903.3680 849.4264 914.6573
## Sep 2018       875.3210 850.6907 899.9513 837.6522 912.9898
## Oct 2018       869.2970 841.3616 897.2324 826.5735 912.0205
## Nov 2018       863.9973 832.5220 895.4726 815.8600 912.1346
## Dec 2018       859.1129 823.6968 894.5290 804.9486 913.2771
## Jan 2019       854.1979 814.4007 893.9952 793.3333 915.0626
## Feb 2019       848.8962 804.3659 893.4264 780.7931 916.9993
## Mar 2019       843.0886 793.6285 892.5488 767.4458 918.7314
## Apr 2019       836.9080 782.4587 891.3574 753.6349 920.1811
## May 2019       830.6354 771.2022 890.0685 739.7403 921.5304
## Jun 2019       824.5462 760.1182 888.9742 726.0121 923.0803
## Jul 2019       818.7853 749.2819 888.2886 712.4891 925.0815
## Aug 2019       813.3226 738.5839 888.0613 699.0196 927.6255
## Sep 2019       807.9966 727.8103 888.1830 685.3622 930.6310
## Oct 2019       802.6108 716.7582 888.4633 671.3106 933.9109
## Nov 2019       797.0291 705.3270 888.7312 656.7828 937.2754
## Dec 2019       791.2284 693.5475 888.9094 641.8383 940.6186
## Jan 2020       785.2897 681.5465 889.0330 626.6281 943.9513
pred.tr <- predict(fit, n.ahead = 12)
U.tr= pred.tr$pred + 2*pred.tr$se
L.tr= pred.tr$pred - 2*pred.tr$se

#Full plot
ts.plot(sqrt(data_testing_ts), xlim = c(2010, 2019), ylim=c(750, 1500), main = "Transformed Unemployment with Forecasts")
lines(U.tr, col="lightblue", lty="dashed")
lines(L.tr, col="lightblue", lty="dashed")
points(seq(from = 2018, to = 2019, length.out=12), pred.tr$pred, col="orchid")

#Zoomed plot 
ts.plot(sqrt(data_testing_ts), xlim = c(2017, 2019), ylim=c(750, 1500), main = "Figure 5: Transformed Unemployment with Forecasts (Zoomed)")
lines(U.tr, col="lightblue", lty="dashed")
lines(L.tr, col="lightblue", lty="dashed")
points(seq(from = 2018, to = 2019, length.out=12), pred.tr$pred, col="orchid")

#We can clearly see that the transformed data from 2019 falls within the prediction interval even if the forecasts are not exactly on the data.
#Original data plot with 12 forecasts
#Full plot
pred.orig <- pred.tr$pred ^ 2
U = U.tr ^ 2
L = L.tr ^ 2
ts.plot(data_testing_ts, xlim = c(2010, 2019), ylim=c(100000, 2500000), main = "Actual Unemployment with Forecasts")
lines(U, col="lightblue", lty="dashed")
lines(L, col="lightblue", lty="dashed")
points(seq(from = 2018, to = 2019, length.out=12), pred.orig, col="orchid")

#Zoomed plot
pred.orig <- pred.tr$pred ^ 2
U = U.tr ^ 2
L = L.tr ^ 2
ts.plot(data_testing_ts, xlim = c(2017, 2019), ylim=c(600000, 1200000), main = "Figure 6: Actual Unemployment with Forecasts (Zoomed)")
lines(U, col="lightblue", lty="dashed")
lines(L, col="lightblue", lty="dashed")
points(seq(from = 2018, to = 2019, length.out=12), pred.orig, col="orchid")

#We can clearly see that the actual data from 2019 falls within the prediction interval even if the forecasts are not exactly on the data.